Thalamocortical Connectivity: Atlas Characteristics and Tract Anatomy

Thalamocortical Autotrack Atlas Anatomical Characteristics

Glasser region names/labels

glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv") #mesulam assignments for glasser parcels
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)

Spatial permutations of the glasser parcellation for spin testing

#https://github.com/frantisekvasa/rotate_parcellation/tree/master, cloned 8/5/2023
source("/cbica/projects/thalamocortical_development/software/rotate_parcellation/R/rotate.parcellation.R")
source("/cbica/projects/thalamocortical_development/software/rotate_parcellation/R/perm.sphere.p.R")
glasser.coords <- read.table("/cbica/projects/thalamocortical_development/software/rotate_parcellation/sphere_HCP.txt", header=F) #The spherical centroids of the multimodal HCP parcellation on the freesurfer sphere; the 360 regions are originally ordered here as left hemisphere -> right hemisphere. The glasser region list and the S-A axis are ordered right hemisphere -> left hemisphere however so let's update the ordering here for spin tests to ensure correspondence
glasser.coords <- rbind(glasser.coords[181:360,], glasser.coords[1:180,])

perm.id.full <- rotate.parcellation(coord.l = as.matrix(glasser.coords[181:360,]), coord.r = as.matrix(glasser.coords[1:180,]), nrot = 10000, method = "vasa") #rotate the glasser parcellation 10,000 times on the freesurfer sphere to generate spatial nulls for spin-based permutation significance testing 
saveRDS(perm.id.full, "/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds")
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins for spin testing (spherical permutation based p-values). Parcel order is right hemisphere -> left hemisphere 
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Glasser parcel anatomy: surface area and sulcal depth measures

glasserparcel.anatomy <- read.csv("/cbica/projects/thalamocortical_development/Templates/fsaverage/stats/glasserparcel_anatomy.csv") %>% select(StructName, Area_mm2, Mean) %>% filter(StructName != "???") %>% set_names("orig_parcelname", "surface_area", "sulcal_depth") #parcel name, total surface area, and mean sulcal depth for 360 glasser regions, calculated in /parcel_anatomy/glasser_parcel_anatomy bash and python scripts
glasserparcel.anatomy <- merge(glasserparcel.anatomy, glasser.assignments, by = "orig_parcelname", sort = F)

Thalamocortical connection atlas tract list

tc.autotrack.tracts <- read.csv("/cbica/projects/thalamocortical_development/software/thalamocortical_autotrack_template/dsi-studio/atlas/ICBM152_adult/ICBM152_adult.tt.gz.txt", header = F) #238 connections
colnames(tc.autotrack.tracts) <- c("tract")
tc.autotrack.tracts$tract <- gsub("-", "_", tc.autotrack.tracts$tract) #no dasher -, no prancer, no vixen

Spin-based enrichment testing function

#A function to test whether an anatomical measure of interest is significant larger or smaller in cortical regions without thalamic connections included in the atlas, based on a spin-based null distribution of the measure
source("/cbica/projects/thalamocortical_development/code/thalamocortical_development/results/tract_anatomy/tractometryatlas_enrichment.R")

Cortical Surface Coverage

Percent of the cortical surface reached by thalamocortical atlas connections

#total area of the cortical surface (sum of the area of all glasser parcels)
total.surfacearea <- sum(glasserparcel.anatomy$surface_area)

#area of the cortical surface with reconstructed thalamus connections
tract.surfacearea <- sum(glasserparcel.anatomy[glasserparcel.anatomy$tract %in% tc.autotrack.tracts$tract,]$surface_area)

#percent of surface area with thalamus-cortex structural connections
round((tract.surfacearea/total.surfacearea)*100,2)
## [1] 76.82

Anatomical Correlates of Reconstructed Connections

Surface area of parcels with versus without reconstructed thalamocortical connections

glasserparcel.anatomy$reconstructed <- ifelse(glasserparcel.anatomy$tract %in% tc.autotrack.tracts$tract, "Present", "Absent")
glasserparcel.anatomy$reconstructed <- factor(glasserparcel.anatomy$reconstructed, ordered = T, levels = c("Present", "Absent"))
glasserparcel.anatomy <- glasserparcel.anatomy %>% mutate(surface_area_z = scale(surface_area, center = T, scale = T)) #z-score surface area
ggplot(glasserparcel.anatomy, aes(x = reconstructed, y = surface_area_z)) + 
  geom_point(aes(color = reconstructed, fill = reconstructed), position = position_jitter(width = .4), size = 0.07) + 
  geom_boxplot(aes(color = reconstructed), fill = alpha("white", 0.75), size = 0.15, outlier.shape = NA) + 
  theme_minimal() +
  scale_color_manual(values = c("#48457F", "#FCAB6A")) +
  labs(x="\n", y="surface area\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.major = element_line(linewidth = 0.2),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1/Atlasconnections_surfacearea.pdf", device = "pdf", dpi = 500, width = 1.3 , height = 1.2)
atlas.enrichment(measure = "surface_area", enrichment_type = "smaller")
## [1] 0

Sulcal depth of parcels with versus without reconstructed thalamocortical connections

#calculate percentage of maximal sulcal depth
#in freesurfer, gyri have negative values and sulci have positive values as they reflect distance from a mid-surface of 0. Let's turn this into percent of maximum sulcal depth
total_sulcalrange <- max(glasserparcel.anatomy$sulcal_depth) - min(glasserparcel.anatomy$sulcal_depth) #span of sulcal depth data
distance_from_gyralcrown <- glasserparcel.anatomy$sulcal_depth - (min(glasserparcel.anatomy$sulcal_depth)) #relative distance from gyral crown
glasserparcel.anatomy$percent_sulcaldepth <- (distance_from_gyralcrown/total_sulcalrange)*100 #percent of total depth
ggplot(glasserparcel.anatomy, aes(x = reconstructed, y = percent_sulcaldepth)) + 
  geom_point(aes(color = reconstructed, fill = reconstructed), position = position_jitter(width = .4), size = 0.07) + 
  geom_boxplot(aes(color = reconstructed), fill = alpha("white", 0.75), size = 0.25, outlier.shape = NA) + 
  theme_minimal() +
  scale_color_manual(values = c("#48457F", "#FCAB6A")) +
  labs(x="\n", y="Percentage sulcal depth\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.major = element_line(linewidth = 0.2),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1/Atlasconnections_sulcaldepth.pdf", device = "pdf", dpi = 500, width = 1.3 , height = 1.2)
atlas.enrichment(measure = "percent_sulcaldepth", enrichment_type = "larger")
## [1] 0.0348

Thalamic Nuclei-specific Cortical Connection Zones

#Number of voxels occupied by each thalamocortical connection in each thalamic nucleus, generated by /thalamocortical_structuralconnectivity/template/thalamocortical_connection_nuclei scripts and based on the Morel-guided MR atlas from Saranathan et al., 2021, Scientific Data
nucleus.connection.profiles <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/template/HCP-MMP_thalamicnuclei_terminationzones.csv") %>% select(tract, label2, label4, label5, label6, label7, label8, label9, label10, label11, label12) #labels 1 and 3 are not used *"unlabeled") in this atlas
colnames(nucleus.connection.profiles) <- c("tract", "AV", "VA", "VLa", "VLp", "VPL", "Pul", "LGN", "MGN", "CM", "MD") #label-nucleus assignments from /cbica/projects/thalamocortical_development/Maps/thalamicnuclei_atlas_saranathan/CustomAtlas.ctbl
connection.list <- gsub("-", "_", nucleus.connection.profiles$tract)
#A function to plot cortical regions that show strongest connectivity to selected thalamic nuclei
nucleus.connection.plot <- function(nuclei.names, n.regions, cortex.color){
    nucleus.connection.profiles %>% select(all_of(nuclei.names)) %>% #select thalamic nuclei
    rowSums %>% as.data.frame() %>% set_names("connectivity.strength") %>% #sum connectivity counts 
    mutate(tract = connection.list) %>% merge(., glasser.regions, by = "tract") %>% 
    slice_max(order_by = connectivity.strength, n = n.regions, with_ties = F) %>% #select top connected cortical regions
    mutate(top.connected = "yes") %>% 
    select(label, top.connected) %>% mutate(cortex = "cortex") %>% 
    rbind(., data.frame(label = "rh_???", top.connected = "no", cortex = "medialwall")) %>%
    rbind(., data.frame(label = "lh_???", top.connected = "no", cortex = "medialwall")) %>%
    filter(grepl("lh_", label)) %>% 
    ggplot(.) + 
      geom_brain(atlas = glasser, hemi = "left", mapping = aes(fill = top.connected, colour=I("black"), size=I(.05))) +
      theme_void() + 
      scale_fill_manual(values = c("white", cortex.color), na.value = "white") + 
      new_scale_fill() + 
      geom_brain(atlas = glasser, hemi = "left", mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
      scale_fill_manual(values = c(alpha("white", 0), "grey"), na.value = "white") +
      theme(legend.position = "none")
}

VLa and VLp

nucleus.connection.plot(nuclei.names = c("VLa", "VLp"), n.regions = 30, cortex.color = "#ffc148")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/VLaVLp_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

VPL

nucleus.connection.plot(nuclei.names = c("VPL"), n.regions = 20, cortex.color = "#faaa6b")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/VPL_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

Pulvinar, LGN

nucleus.connection.plot(nuclei.names = c("Pul", "LGN"), n.regions = 30, cortex.color = "#e4948b")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/PulLGN_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

CM

nucleus.connection.plot(nuclei.names = c("CM"), n.regions = 25, cortex.color = "#b65099")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/CM_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

MD

nucleus.connection.plot(nuclei.names = c("MD"), n.regions = 55, cortex.color = "#7f3685")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/MD_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

AV

nucleus.connection.plot(nuclei.names = c("AV"), n.regions = 50, cortex.color = "#47378f")
## merging atlas and data by 'label'
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/AV_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)
## merging atlas and data by 'label'
## merging atlas and data by 'label'

Thalamocortical Connectivity Cortical and Thalamic Gradients

S-A axis ranks

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by="orig_parcelname", sort = F)

Dataset-specific final tract lists

#Tract lists generated by /tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv")
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv")

Function to calculate tract-level means per dataset

tract_average_measures <- function(measure, atlas, dataset){
  input.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset)) %>% select(contains("autotrack")) #nrow subjects by ncol tracts for mean calculation
  tract.means <- colMeans(input.df, na.rm = TRUE) %>% as.data.frame %>% set_names(sprintf("mean_%s", measure)) %>% mutate(tract = row.names(.))
  tract.means <- merge(tract.means, S.A.axis, by = "tract")
  return(tract.means)
}

Thalamocortical Connection Fractional Anisotropy

FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #FA for included tracts and participants

FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv")

Calculate mean FA for each connection in PNC and HCPD datasets

tract.meanFA.pnc <- tract_average_measures(measure = "FA", atlas = "glasser", dataset = "pnc")
tract.meanFA.pnc <- tract.meanFA.pnc[tract.meanFA.pnc$tract %in% tractlist.pnc$tract,] #connection-level mean FA for final PNC tract list
write.csv(tract.meanFA.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_FA_primary.csv", row.names = F, quote = F)

tract.meanFA.hcpd <- tract_average_measures(measure = "FA", atlas = "glasser", dataset = "hcpd")
tract.meanFA.hcpd <- tract.meanFA.hcpd[tract.meanFA.hcpd$tract %in% tractlist.hcpd$tract,] #connection-level mean FA for final HCPD tract list
write.csv(tract.meanFA.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_FA_primary.csv", row.names = F, quote = F)

Correlation of connection-level average FA between datasets

Pearson’s r

tract.meanFA.merged <- merge((tract.meanFA.pnc %>% select(tract, mean_FA)), (tract.meanFA.hcpd %>% select(tract, mean_FA)), by = "tract") %>% set_names("tract", "mean_FA.pnc", "mean_FA.hcpd")

cor.test(tract.meanFA.merged$mean_FA.pnc, tract.meanFA.merged$mean_FA.hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  tract.meanFA.merged$mean_FA.pnc and tract.meanFA.merged$mean_FA.hcpd
## t = 61.101, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9635117 0.9783814
## sample estimates:
##       cor 
## 0.9719002

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, tract.meanFA.merged, by = "tract") #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$mean_FA.pnc, y = spin.data$mean_FA.hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0

Correlation plot

ggplot(tract.meanFA.merged, aes(x = mean_FA.pnc, y = mean_FA.hcpd)) +
  geom_point(color = c("#571766"), size = 0.5) +
  geom_smooth(method = "lm", color = c("goldenrod1"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(0.2, 0.525)) +
  scale_y_continuous(limits = c(0.2, 0.525)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/FA_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.4)

Thalamocortical connection FA organization along the S-A axis

Spearman’s correlations

PNC

cor.test(tract.meanFA.pnc$mean_FA, tract.meanFA.pnc$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tract.meanFA.pnc$mean_FA and tract.meanFA.pnc$SA.axis
## S = 2633482, p-value = 4.941e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4248717

HCPD

cor.test(tract.meanFA.hcpd$mean_FA, tract.meanFA.hcpd$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tract.meanFA.hcpd$mean_FA and tract.meanFA.hcpd$SA.axis
## S = 2877906, p-value = 4.757e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4192293

Spatial permutation (spin) based p-value

PNC

spin.data <- left_join(spin.parcels, tract.meanFA.pnc, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_FA, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.045

HCPD

spin.data <- left_join(spin.parcels, tract.meanFA.hcpd, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_FA, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.06135

Thalamocortical FA by S-A axis plot

PNC

ggplot(tract.meanFA.pnc, aes(x = SA.axis, y = mean_FA)) +
  geom_point(aes(color = SA.axis), size = 0.5) +
  geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) + 
  theme_classic() +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
  labs(x="\nSensorimotor-association axis rank", y="Connection FA\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/FA_SArank_PNC.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)

HCPD

ggplot(tract.meanFA.hcpd, aes(x = SA.axis, y = mean_FA)) +
  geom_point(aes(color = SA.axis), size = 0.5) +
  geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) + 
  theme_classic() +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
  labs(x="\nSensorimotor-association axis rank", y="Connection FA\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/FA_SArank_HCPD.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)

Thalamocortical Connection Core-Matrix Gradient (CPt) Values

Read in connection-specific core-matrix gradient values (based on the connection’s thalamic endpoints and the thalamic calbindin-parvalbumin CPt gradient) and get data for just final study sample and included tracts

extract_measure <- function(metric.df, participants.df, measure){ 
  measure.df <- metric.df %>% select(rbcid, tract, all_of(measure)) #extract measure of interest from long df
  measure.df <- measure.df %>% arrange(tract) #tracts in alphabetical order 
  measure.df <- measure.df %>% pivot_wider(names_from = "tract", values_from = all_of(measure)) #long to wide format, all participants and tracts
  measure.df <- measure.df %>% arrange(rbcid)  #rbcids in numerical order
  colnames(measure.df) <- gsub("-", "_", colnames(measure.df))
  measure.df <- measure.df[measure.df$rbcid %in% participants.df$rbcid,] #only final sample of participants now
  return(measure.df)
}

PNC

#Create connection CPt value spreadsheet for final sample of PNC participants
CPt.glasser.pnc.long <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_thalamicconnection_CPt.csv") #thalamic connection CPt gradient values, generated by /tract_measures/CPtgradient_values/lthalamocortical_calculate_CPtvalues.sh and /tract_measures/CPtgradient_values/thalamocortical_CPtvalues.py
participants.pnc <- read.csv("/cbica/projects/thalamocortical_development/sample_info/PNC/PNC_thalamocortical_finalsample_N1145.csv") #PNC final participant list
CPt.glasser.pnc <- extract_measure(metric.df = CPt.glasser.pnc.long, participants.df = participants.pnc, measure = "CPt")
#Don't analyze CPt values for connections with <5 streamlines 
streamlinecount.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_streamlinecount_finalsample.csv")
identical(streamlinecount.glasser.pnc$rbcid, CPt.glasser.pnc$rbcid) 
## [1] TRUE
streamlinecount.glasser.pnc <- streamlinecount.glasser.pnc %>% select(contains("autotrack")) #1145 subjects by 238 tracts
CPt.glasser.pnc <- CPt.glasser.pnc %>% select(contains("autotrack")) #1145 subjects by 238 tracts
identical(names(streamlinecount.glasser.pnc), names(CPt.glasser.pnc))
## [1] TRUE
CPt.glasser.pnc <- CPt.glasser.pnc * ifelse(streamlinecount.glasser.pnc >= 5, 1, NA)

HCPD

#Create connection CPt value spreadsheet for final sample of HCPD participants
CPt.glasser.hcpd.long <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_thalamicconnection_CPt.csv") 
participants.hcpd <- read.csv("/cbica/projects/thalamocortical_development/sample_info/HCPD/HCPD_thalamocortical_finalsample_N572.csv") #HCPD final participant list
CPt.glasser.hcpd <- extract_measure(metric.df = CPt.glasser.hcpd.long, participants.df = participants.hcpd, measure = "CPt")
#Don't analyze CPt values for connections with <5 streamlines 
streamlinecount.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_streamlinecount_finalsample.csv")
identical(streamlinecount.glasser.hcpd$rbcid, CPt.glasser.hcpd$rbcid)
## [1] TRUE
streamlinecount.glasser.hcpd <- streamlinecount.glasser.hcpd %>% select(contains("autotrack")) #572 subjects by 238 tracts
CPt.glasser.hcpd <- CPt.glasser.hcpd %>% select(contains("autotrack")) #572 subjects by 238 tracts
identical(names(streamlinecount.glasser.hcpd), names(CPt.glasser.hcpd)) 
## [1] TRUE
CPt.glasser.hcpd <- CPt.glasser.hcpd * ifelse(streamlinecount.glasser.hcpd >= 5, 1, NA)

Calculate thalamic connection area CPt value for each tract in PNC and HCPD datasets

tract.meanCPt.pnc <- tract_average_measures(measure = "CPt", atlas = "glasser", dataset = "pnc")
tract.meanCPt.pnc <- tract.meanCPt.pnc[tract.meanCPt.pnc$tract %in% tractlist.pnc$tract,] 
write.csv(tract.meanCPt.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv", row.names = F, quote = F)

tract.meanCPt.hcpd <- tract_average_measures(measure = "CPt", atlas = "glasser", dataset = "hcpd")
tract.meanCPt.hcpd <- tract.meanCPt.hcpd[tract.meanCPt.hcpd$tract %in% tractlist.hcpd$tract,] 
write.csv(tract.meanCPt.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv", row.names = F, quote = F)

Correlation of connection-level average CPt between datasets

Pearson’s r

tract.meanCPt.merged <- merge((tract.meanCPt.pnc %>% select(tract, mean_CPt)), (tract.meanCPt.hcpd %>% select(tract, mean_CPt)), by = "tract") %>% set_names("tract", "mean_CPt.pnc", "mean_CPt.hcpd")

cor.test(tract.meanCPt.merged$mean_CPt.pnc, tract.meanCPt.merged$mean_CPt.hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  tract.meanCPt.merged$mean_CPt.pnc and tract.meanCPt.merged$mean_CPt.hcpd
## t = 93.485, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9839914 0.9905554
## sample estimates:
##       cor 
## 0.9877012

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, tract.meanCPt.merged, by = "tract") 
perm.sphere.p(x = spin.data$mean_CPt.pnc, y = spin.data$mean_CPt.hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0

Correlation plot

ggplot(tract.meanCPt.merged, aes(x = mean_CPt.pnc, y = mean_CPt.hcpd)) +
  geom_point(color = c("#571766"), size = 0.5) +
  geom_smooth(method = "lm", color = c("goldenrod1"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/CPt_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.4)

Thalamocortical connection CPt organization along the S-A axis

Spearman’s correlations

PNC

cor.test(tract.meanCPt.pnc$mean_CPt, tract.meanCPt.pnc$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tract.meanCPt.pnc$mean_CPt and tract.meanCPt.pnc$SA.axis
## S = 781550, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5771346

HCPD

cor.test(tract.meanCPt.hcpd$mean_CPt, tract.meanCPt.hcpd$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tract.meanCPt.hcpd$mean_CPt and tract.meanCPt.hcpd$SA.axis
## S = 1015854, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4990352

Spatial permutation (spin) based p-value

PNC

spin.data <- left_join(spin.parcels, tract.meanCPt.pnc, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_CPt, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.00835

HCPD

spin.data <- left_join(spin.parcels, tract.meanCPt.hcpd, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_CPt, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.04045

Thalamocortical FA by S-A axis plot

PNC

ggplot(tract.meanCPt.pnc, aes(x = SA.axis, y = mean_CPt)) +
  geom_point(aes(color = SA.axis), size = 0.5) +
  geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) + 
  theme_classic() +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
  labs(x="\nSensorimotor-association axis rank", y="Connection CPt\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/CPt_SArank_PNC.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)

HCPD

ggplot(tract.meanCPt.hcpd, aes(x = SA.axis, y = mean_CPt)) +
  geom_point(aes(color = SA.axis), size = 0.5) +
  geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) + 
  theme_classic() +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
  labs(x="\nSensorimotor-association axis rank", y="Connection CPt\n") +
  theme(legend.position = "none") +
  theme(
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/CPt_SArank_HCPD.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)